      subroutine hplset(xxlp,dur,dynam,plamp,seed,sr,q)                  
c...............need to replace rand unit              
c.................compile with order -lF77 -lm    
      dimension q(1)                                                    
      complex j                                                         
c------------------decay stretch & shortening factors                   
c	write(6,999)xlp,dur,dynam,plamp,seed,sr
c999   format(6f10.4)
      xlp=xxlp
      freq=1./xlp                                                       
      pi=3.141592654
      w=2.*pi*freq                                                      
c.....tau => decay time required                                        
      tau=dur/6.91                                                      
      n1=aint(sr/freq-.5)                                               
21    ga1=exp(-(n1+.5)/(tau*sr))                                        
      x=1-ga1**2.                                                       
      b=2.*(1.-cos(w/sr))                                               
      s2=b**2.-4.*b*x                                                   
      if(s2.lt.0.) goto 11                                              
c------------------------                                               
      s1=sqrt(s2)/(2*b)                                                 
      s1=.5+s1                                                          
c.....rho => loss factor                                                
      rho=1.0                                                           
      if(s1.lt.1) goto 10                                               
c------------------------                                               
11    s1=.5                                                             
      rho=exp(-1./(tau*freq))/abs(cos(w/(2*sr)))                        
c------------------------                                               
10    c=1.-s1                                                           
      c1=s1                                                             
c----------------- tuning filter coefficient                            
      xlp=sr/freq                                                       
c.....pa => phase delay due to the 2 point average filter               
      pa=(-1)*sr/w*atan(-c*sin(w/sr)/((1.-c)+c*cos(w/sr)))              
      pa=1.-pa                                                          
c.....n => buffer length                                                
      n=aint(xlp-pa)                                                    
c.....ga=sqrt((c1**2)+(c**2)+2.*c*c1*cos(w/sr))                         
c.....tau=-(n+pa)/(sr*alog(ga))                                         
      if(n.eq.n1) goto 20                                               
      n1=n                                                              
      goto 21                                                           
c.....pc => delay required from the allpass filter                      
20    pc=(sr/freq)-n-pa                                                 
c.....cc => coefficient of the allpass filter                           
      cc=sin((w/sr-w/sr*pc)/2.)/sin((w/sr+w/sr*pc)/2.)                  
c........                                                               
c----------------- dynamic filter coefficient                           
      fl=20.                                                            
      fu=sr/2.                                                          
      fm=sqrt(fl*fu)                                                    
      rl=exp(-pi*dynam/sr)                                              
      j=cmplx(0.,-1.)                                                   
      gl=(1.-rl)/cabs(1.-(rl*cexp(j*2*pi*fm/sr)))                       
      gl2=gl**2                                                         
      w2=w/2                                                            
      r1=(1.-gl2*cos(w/sr))/(1-gl2)                                     
      r2=2*gl*sin(w2/sr)*sqrt(1.-gl2*(cos(w2/sr)**2))/(1.-gl2)          
      r=r1+r2                                                           
      if(r.ge.1) r=r1-r2                                                
c-------------------                                                    
      q(2)=n+20                                                         
      len=q(2)                                                          
      do 1 m=20,len                                                     
        seed=frand(seed)
	x=seed*2-1.
        q(m)=plamp                                                      
        if(x.lt.0.) q(m)=-plamp                                         
1     continue                                                          
c--------------------                                                   
	q(2)=q(2)-1
c.......reset for C code which starts from 0
      q(1)=q(2)                                                         
      q(3)=c                                                            
      q(4)=1.-c                                                         
      q(9)=r                                                            
      q(10)=rho                                                         
      q(11)=cc                                                          
      q(5)=0.                                                           
      q(6)=0.                                                           
      q(7)=0.                                                           
      q(8)=0.                                                           
      return                                                            
      end                                                               
      function frand(x)
      n=x*1048576.
      frand=float(mod(1061*n+221589,1048576))/1048576.
      return
      end
